RangeShiftR Workshop
Anne-Kathleen Malchow, Greta Bocedi, Justin Travis, Damaris Zurell
October 21, 2019
In this tutorial, we will revisit the examples provided in the original RangeShifter publication (Bocedi et al. 2014). This will introduce you to the main features of the new R package RangeShiftR. We will use the same parameters as in the original publication and the accompanying tutorial based on the Windows GUI. Also, we provide codes for plotting figures similar to those provided in the Windows GUI and in (Bocedi et al. 2014).
1 Installing RangeShiftR
The RangeshiftR package is not yet released and we thus need to install it by hand.
1.1 Install and load required packages
RangeshiftR requires R 3.6.0 or higher and Rcpp 1.0.1 or higher. Devtools is recommended for an easy installation. Install these packages if you don’t have them already.
library(Rcpp)
library(devtools)## Loading required package: usethis1.2 Compiler toolchain
In order to build the package from source, we need to set up a compiler toolchain (i.e. install a suitable C++ compiler and make it available to R).
on Windows: the Rtools package on CRAN contains all the required tools. Find installing instructions - for either manual installation or from within R (further down) - here: https://thecoatlessprofessor.com/programming/installing-rtools-for-compiled-code-via-rcpp/
Pay extra attention to steps 10/11. In case you have to add the path variable manually, try this (Windows 10):
Go to Control panel -> System and Security -> System -> Advanced system settings -> Environmental variables -> in the list of system variables, edit ‘Path’ -> add the 2 variables ‘C:\Rtools\bin’, ‘C:\Rtools\mingw_32\bin’on Linux: Install the standard development packages (e.g. on Ubuntu and Debian: r-base-dev).
on MacOS: Install the Apple Developer Tools Xcode and the C++ compiler clang7. Find installing instructions here (for RangeShifter, the Fortran compiler gfortran is not needed, so you may leave this step out): https://thecoatlessprofessor.com/programming/r-compiler-tools-for-rcpp-on-macos/
You can check if everything works properly by trying this in your R terminal:
Rcpp::evalCpp("2+2")## [1] 4If the output is not a number (i.e. 4), then there is an issue with the toolchain.
1.3 Install the package
Navigate to package directory “RangeshiftR” and run the installation.
devtools::install()After successful installation, you can load the package:
library(RangeshiftR)2 Examples & Tutorials
In this section we will go step by step through the first two examples presented in the companion paper (Bocedi et al. 2014). These will cover some of the main features of RangeShifter, and help in becoming familiar with the software. We will use the same parameters that were used in the paper. Nevertheless, we encourage you do to more experimenting and trying different parameters and combinations of options for getting to know the program.
2.1 Exercise 1:
Range expansion, long-distance dispersal and environmental stochasticity
This is an example of how RangeShifter can be used at national scale for modelling species range dynamics. Here we model a hypothetical grassland species distributed initially in the South-West of England, and assume that from the start of the simulation the species is free to expand its range. This could be the case for alien species that naturally start to expand after having gone through an establishment phase, alien or native species that have been released from natural enemies or competitors, or species for which a previously prohibiting climate has become suitable. We assume that we have data about the current species distribution and use it as a starting point. The objective is to investigate how different assumptions about the dispersal ability of the species and about temporal environmental stochasticity can affect the modelled range expansion. We start with the basic setting of a single dispersal kernel and no environmental stochasticity.
2.1.1 Create a RS directory
The directory in which we run RangeshifteR needs to have a certain folder structure. It should contain 3 sub-folders named ‘Inputs’, ‘Outputs’ and ‘Output_Maps’. We can create them from R:
# relative path from working directory:
dirpath = "Tutorials/Exercise1/"
# dir.create(paste0(dirpath,"Inputs"), showWarnings = TRUE)
# dir.create(paste0(dirpath,"Outputs"), showWarnings = TRUE)
# dir.create(paste0(dirpath,"Output_Maps"), showWarnings = TRUE)Copy the map files provided for exercise 1 into the ‘Inputs’ folder.
2.1.2 Landscape parameters
We use a land-cover map of Great Britain at 1km resolution. Six dominant aggregated habitat types were derived from LandCover Map 2007 (Morton et al. 2011). The map, UKmap_1km.txt, is a ASCII in the standard text format, where each cell holds the code of its dominant habitat type. For simplicity, the codes are set as sequential numbers:
- 1 = woodland (broadleaved and conifer)
- 2 = arable
- 3 = improved grassland
- 4 = semi-natural grassland (acid, neutral and calcareous grassland)
- 5 = heath and bog
- 6 = other (urban, water & coastal habitats)
Furthermore, we are provided with a map that defines the (initial) species distribution, named Species_Distribution_10km.txt. It is given at a larger resolution of 10km.
We can plot the landscape map and overlay it with the initial species distribution:
library(raster)
library(RColorBrewer)
library(rasterVis)
library(latticeExtra)
library(viridis)
library(grid)
library(gridExtra)
UKmap <- raster(paste0(dirpath, "Inputs/UKmap_1km.txt"))
SpDist <- raster(paste0(dirpath, "Inputs/Species_Distribution_10km.txt"))
values(SpDist)[values(SpDist) < 1] <- NA
# Plot land cover map and highlight cells with initial species distribution - option 1:
plot(UKmap, col=brewer.pal(n = 6, name = "Spectral"), axes=F)
plot(rasterToPolygons(SpDist, dissolve=F), add=T)
# Plot land cover map and highlight cells with initial species distribution - option 2 with categorical legend:
UKmap.f <- as.factor(UKmap)
# add the land cover classes to the raster attribute table (rat)
rat <- levels(UKmap.f)[[1]]
rat[["landcover"]] <- c("woodland", "arable", "improved grassland", "semi-natural grassland", "heath and bog", "other")
levels(UKmap.f) <- rat
custom.pal <- c("#1A9850", "#91CF60", "#D9EF8B", "#FEE08B", "#D8B365", "#777777")
levelplot(UKmap.f, margin=F, scales=list(draw=FALSE), col.regions=custom.pal) +
layer(sp.polygons(rasterToPolygons(SpDist, dissolve=F), fill=NA, col='black'))
In order to use the habitat and distribution maps in RangeShiftR, we need to set up a landscape parameter object. It takes the path to the map files along with some other parameters, like their resolutions (given in meters) and the number of habitats. Additionally, for each habitat type we have to specify the carrying capacity for our target species, given in individuals per hectare. In this example, we assume that the species can reproduce only in semi-natural grassland, which has the code 4. Set the carrying capacity for this habitat equal to 5 inds/ha and all others to 0, by defining a vector (carrycap) that contains these numbers in order of increasing habitat codes:
# carrying capacitíes and landscape parameter object
carrycap <- c(0, 0, 0, 5, 0, 0)
land <- ImportedLandscape(LandscapeFile = "UKmap_1km.txt",
Resolution = 1000,
Nhabitats = 6,
K = carrycap,
SpDistFile = "Species_Distribution_10km.txt",
SpDistResolution = 10000)2.1.3 Species parameters
Next, we specify the species parameters for demography and dispersal kernel. In this first example, we choose the most simple model options for both: the population dynamics are described by an only-female model with non-overlapping generations. It has one compulsory argument, the maximum growth rate Rmax, that we set to 1.5.
demo <- Demography(Rmax = 1.5)The dispersal module comprises three sub-modules for Emigration, Transport and Settlement. We assume a constant emigration probability of 0.1. The transfer phase is modelled with a dispersal kernel, whose mean distance we set to 2,000m. For the Settlement, we keep all default options, which means that if an individual arrives into an unsuitable cell it dies.
disp <- Dispersal(Emigration = Emigration(EmigProb = 0.1),
Transfer = DispersalKernel(Distances = 2000),
Settlement = Settlement() )
# Alternative notation:
# disp <- Dispersal() + Emigration(EmigProb = 0.1) + DispersalKernel(Distances = 2000) + Settlement()2.1.4 Initialisation parameters
In order to control the initial distribution of individuals in the landscape at year 0, we set initialisation rules. We want to initialise all the habitat cells that are present inside the 10km x 10km cells of the loaded species distribution map; and we want each of those cells to be initialised at its carrying capacity. These options are encoded in numeric arguments, for a list of possible settings see the documentation:
?Initialise## No documentation for 'Initialise' in specified packages and libraries:
## you could try '??Initialise'init <- Initialise(InitType = 1, # = initialisation from a loaded species distribution map
SpType = 0, # = all suitable cells within all distribution presence cells
InitDens = 0) # = at carrying capacity2.1.5 Simulation parameters
Lastly, we set some basic simulation parameters, i.e. the simulation number, number of replicates and number of years to be simulated. Furthermore, we specify what file output will be generated. In this example, we choose to output the population, occupancy and range data for every 5 years.
sim_0 <- Simulation(Simulation = 0,
Replicates = 20,
Years = 100,
OutIntPop = 5,
OutIntOcc = 5,
OutIntRange = 5)2.1.6 Parameter master
All the settings we have made so far are contained in their respective modules. They have to be combined in a parameter master object, which is needed to run the simulation.
s <- RSsim(land = land, demog = demo, dispersal = disp, simul = sim_0, init = init)
# Alternative notation:
# s <- RSsim() + land + demo + disp + sim_0 + initBefore we run the simulation, let’s get an overview of the settings we have made by simply typing
s## Simulation # 0
## -----------------
## Replicates = 20
## Years = 100
## Absorbing = FALSE
## File Outputs:
## Range, every 5 years
## Occupancy, every 5 years
## Populations, every 5 years, starting year 0
##
## Landscape imported from file:
## UKmap_1km.txt
## with 6 unique integer habitat code(s)
## Carrying capacities K = 0 0 0 5 0 0 .
## Resolution : 1000
## Initial Species Distribution imported from file:
## Species_Distribution_10km.txt
## Resolution : 10000
##
## Demography:
## Unstructured population:
## Rmax : 1.5
## bc : 1
## Reproduction Type : 0 (female only)
##
## Dispersal:
## Emigration:
## Emigration probabilities:
## [,1]
## [1,] 0.1
##
## Transfer:
## Dispersal Kernel
## Dispersal kernel traits:
## [,1]
## [1,] 2000
## Constant mortality prob = 0
##
## Settlement:
## Settlement conditions:
## [,1]
## [1,] 0
## FindMate = FALSE
##
## Initialisation:
## InitType = 1 : Initialisation from loaded species distribution map
## all presence cells/patches.
## InitDens = 0 : At carrying capacity KAlso, the RSsim() function (see next section) will check the validity of all settings and give error messages in case something doesn’t match.
2.1.7 Run the simulation
Once all parameters are set and a parameter master object has been defined, we can run the simulations in the specified RS directory.
RunRS(s, dirpath)2.1.8 Plot results
For convenience and for compliance with the RangeShifter Windows GUI, the RangeShiftR package contains a few plotting functions to inspect the population trajectories of the replicate simulations.
range_df <- readRange(s, dirpath)
par(mfrow=c(1,2))
plotAbundance(range_df)
plotOccupancy(range_df)
# Or plot standard deviation:
par(mfrow=c(1,2))
plotAbundance(range_df, rep=F, sd=T)
plotOccupancy(range_df, rep=F, sd=T)
The population output file contains the local abundance for all cells, timesteps and replicates. We can convert this dataframe into raster maps and plot the spatial abundance distribution.
# read population output file into a dataframe
pop_df <- readPop(s, dirpath)
# Not all years have the same number of cells. For later stacking, we need a common extent. This is a quick & dirty solution:
ext <- c(min(pop_df$x)-500,max(pop_df$x)+500,min(pop_df$y)-500,max(pop_df$y)+500)
# Make stack of different raster layers for each year and for only one repetition (Rep==0):
pop_wide_rep0 <- reshape(subset(pop_df,Rep==0)[,c('Year','x','y','NInd')], timevar='Year', v.names=c('NInd'), idvar=c('x','y'), direction='wide')
r_years_rep0 <- rasterFromXYZ(pop_wide_rep0)
# Overlay with UK mask
r_years_rep0 <- extend(r_years_rep0, UKmap)
values(r_years_rep0)[is.na(values(r_years_rep0))] <- 0
r_years_rep0 <- mask(r_years_rep0, UKmap)
# Map abundance
levelplot(r_years_rep0[['NInd.90']], margin=F, scales=list(draw=FALSE), at=c(0,seq(1,max(pop_df$NInd), length=20)),
col.regions=c('grey',rev(inferno(20)))) +
layer(sp.polygons(rasterToPolygons(SpDist, dissolve=F), fill=NA, col='red'))
We can make similar maps showing the average abundance over all replicate runs. For convenience, we define a small function for this, which we can re-use later. When calling this function, we can choose whether all replicates should be averaged (the default) or whether a specific replicate should be extracted. Also, we can optionally overlay a mask, e.g. the UK map.
stack_pop <- function(pop_df, ext, rep=NULL, mask=NULL){
# This function takes the population data frame output from RangeShiftR and turns this into a raster stack of abundance maps.
# If the ID of the Replicate ("rep") is not provided, it will return the mean abundance over all replicates.
if (!is.null(rep)){
pop_wide <- reshape(subset(pop_df,Rep==rep)[,c('Year','x','y','NInd')], timevar='Year', v.names=c('NInd'), idvar=c('x','y'), direction='wide')
r_years <- rasterFromXYZ(pop_wide)
if (!is.null(mask)){
r_years <- extend(r_years, mask)
values(r_years)[is.na(values(r_years))] <- 0
r_years <- mask(r_years, mask)
}
} else {
pop_wide <- lapply(unique(pop_df$Year),FUN=function(year){reshape(subset(pop_df,Year==year)[,c('Rep','x','y','NInd')], timevar='Rep', v.names=c('NInd'), idvar=c('x','y'), direction='wide')})
r_years <- stack(sapply(pop_wide, FUN=function(i){mean(extend(rasterFromXYZ(i),ext))}))
names(r_years) <- paste0('mean.NInd.',unique(pop_df$Year))
if (!is.null(mask)){
r_years <- extend(r_years, mask)
values(r_years)[is.na(values(r_years))] <- 0
r_years <- mask(r_years, mask)
}
}
return(r_years)
}# Extract maps of single replicate run:
r_years_rep0 <- stack_pop(pop_df, ext, rep=0, mask=UKmap)
names(r_years_rep0)## [1] "NInd.0" "NInd.5" "NInd.10" "NInd.15" "NInd.20" "NInd.25"
## [7] "NInd.30" "NInd.35" "NInd.40" "NInd.45" "NInd.50" "NInd.55"
## [13] "NInd.60" "NInd.65" "NInd.70" "NInd.75" "NInd.80" "NInd.85"
## [19] "NInd.90" "NInd.95" "NInd.100"# Extract maps with averaged abundances over all replicates:
r_years <- stack_pop(pop_df, ext, mask=UKmap)
names(r_years)## [1] "mean.NInd.0" "mean.NInd.5" "mean.NInd.10" "mean.NInd.15"
## [5] "mean.NInd.20" "mean.NInd.25" "mean.NInd.30" "mean.NInd.35"
## [9] "mean.NInd.40" "mean.NInd.45" "mean.NInd.50" "mean.NInd.55"
## [13] "mean.NInd.60" "mean.NInd.65" "mean.NInd.70" "mean.NInd.75"
## [17] "mean.NInd.80" "mean.NInd.85" "mean.NInd.90" "mean.NInd.95"
## [21] "mean.NInd.100"# Map abundance
levelplot(r_years_rep0[['NInd.90']], margin=F, scales=list(draw=FALSE), at=c(0,seq(1,max(pop_df$NInd),length=20)), col.regions=c('grey',rev(inferno(20)))) +
layer(sp.polygons(rasterToPolygons(SpDist, dissolve=F), fill=NA, col='red'))
# Map mean abundance
levelplot(r_years[['mean.NInd.90']], margin=F, scales=list(draw=FALSE), at=c(0,seq(1,max(pop_df$NInd),length=20)), col.regions=c('grey',rev(inferno(20)))) +
layer(sp.polygons(rasterToPolygons(SpDist, dissolve=F), fill=NA, col='red'))
2.1.9 Change simulations
We can easily change certain aspects of the simulations by changing the various modules. Here, we explore the different simulations as shown in Fig. 2 of Bocedi et al. (2014). Firstly, we add long-distance dispersal. Secondly, we add environmental stochasticity without and with temporal autocorrelation.
2.1.9.1 Add rare long-distance dispersal
To explore the effect of rare long-distance dispersal events on the range expansion, we only need to change the dispersal module. More specifically, we define a new transfer sub-module with the option DoubleKernel enabled, to use the double negative exponential kernel (Mixed kernel). This means that in addition to the standard dispersal kernel, an individual has a certain probability of dispersing with a second - a long-distance - dispersal kernel.
trans_long <- DispersalKernel(DoubleKernel = T, Distances = matrix(c(2000, 10000, 0.99), ncol = 3)) #Distances = c(MeanDistance-I, MeanDistance-II, probability of dispersing with Kernel-I)
disp_long <- Dispersal(Emigration = Emigration(EmigProb = 0.1),
Transfer = trans_long,
Settlement = Settlement() )Moreover we change the simulation number to 1 in the corresponding parameter object, in order to avoid overwriting the previous output.
sim_1 <- Simulation(Simulation = 1,
Replicates = 20,
Years = 100,
OutIntPop = 5,
OutIntOcc = 5,
OutIntRange = 5)Finally, the new modules are added to the old parameter master to define a new one, s_long.
s_long <- s + disp_long + sim_1With this, run the new simulation:
We can now compare the result plots. The number of occupied cells as well as the total abundance is much higher in this scenario.
# Get range results
range_df_long <- readRange(s_long, dirpath)
# Plot total abundance and number of occupied cells
par(mfrow=c(1,2))
plotAbundance(range_df_long)
plotOccupancy(range_df_long)
# Map mean abundance
pop_df_long <- readPop(s_long, dirpath)
ext_long <- c(min(pop_df_long$x)-500,max(pop_df_long$x)+500,min(pop_df_long$y)-500,max(pop_df_long$y)+500)
r_years_long <- stack_pop(pop_df_long, ext_long, mask=UKmap)
levelplot(r_years_long[['mean.NInd.90']], margin=F, scales=list(draw=FALSE), at=c(0,seq(1,max(pop_df_long$NInd),length=20)), col.regions=c('grey',rev(inferno(20)))) +
layer(sp.polygons(rasterToPolygons(SpDist, dissolve=F), fill=NA, col='red'))
When comparing the maps of mean abundance between the first and second scenario, these look quite similar. One reason might be that long-distance dispersal is highly stochastic and, thus, high occupancy only occurs for a few but not for all runs. So, let’s try and map only the replicate run with the highest occupancy.
# Which replicate had highest occupancy?
head(subset(range_df_long,Year==95)[order(subset(range_df_long,Year==95)$NOccupCells,decreasing=T),])## Rep Year RepSeason NInds NOccupCells Occup.Suit min_X max_X min_Y
## 188 8 95 0 1061177 3342 0.140627 211000 569000 58000
## 62 2 95 0 1126298 3168 0.133305 211000 569000 45000
## 20 0 95 0 924124 2637 0.110961 211000 569000 45000
## 272 12 95 0 785900 2496 0.105028 211000 577000 58000
## 104 4 95 0 687586 2161 0.090932 211000 577000 51000
## 146 6 95 0 726067 2008 0.084494 211000 569000 45000
## max_Y
## 188 328000
## 62 314000
## 20 289000
## 272 315000
## 104 271000
## 146 273000rep_x = subset(range_df_long,Year==95)$Rep[which.max(subset(range_df_long,Year==95)$NOccupCells)]
r_years_long_1 <- stack_pop(pop_df_long, ext_long, rep=rep_x, mask=UKmap)
names(r_years_long_1)## [1] "NInd.0" "NInd.5" "NInd.10" "NInd.15" "NInd.20" "NInd.25"
## [7] "NInd.30" "NInd.35" "NInd.40" "NInd.45" "NInd.50" "NInd.55"
## [13] "NInd.60" "NInd.65" "NInd.70" "NInd.75" "NInd.80" "NInd.85"
## [19] "NInd.90" "NInd.95" "NInd.100"levelplot(r_years_long_1[['NInd.90']], margin=F, scales=list(draw=FALSE), at=c(0,seq(1,max(pop_df_long$NInd),length=20)), col.regions=c('grey',rev(inferno(20)))) +
layer(sp.polygons(rasterToPolygons(SpDist, dissolve=F), fill=NA, col='red'))
2.1.9.2 Add temporal environmental stochasticity
Lastly, we incorporate temporal environmental stochasticity, which is recognized to be fundamental for both ecological and evolutionary processes and is expected to increase in frequency and become more auto-correlated with climate change.
The options for the environmental stochasticity are set in the Simulation parameters object. For this example we use global (referring to the spatial extent) stochasticity in growth rate. We explore two types of stochasticity: temporally uncorrelated (white noise) and positively correlated (red noise), by setting the autocorrelation coefficient ac:
sim_2 <- Simulation(Simulation = 2,
Replicates = 20,
Years = 100,
EnvStoch = 1, # global environmental stochasticity
EnvStochType = 0, # in growth rate:
std = 0.25, # magnitude of stochastic fluctuations
ac = 0.0, # temporal autocorrelation
minR = 0.5, # minimum growth rate
maxR = 2.5, # maximum growth rate
OutIntPop = 5,
OutIntOcc = 5,
OutIntRange = 5)
s_envstoch_white <- s_long + sim_2sim_3 <- Simulation(Simulation = 3,
Replicates = 20,
Years = 100,
EnvStoch = 1, # global environmental stochasticity
EnvStochType = 0, # in growth rate:
std = 0.25, # magnitude of stochastic fluctuations
ac = 0.7, # temporal autocorrelation
minR = 0.5, # minimum growth rate
maxR = 2.5, # maximum growth rate
OutIntPop = 5,
OutIntOcc = 5,
OutIntRange = 5)
s_envstoch_red <- s_long + sim_3RunRS(s_envstoch_white, dirpath)
RunRS(s_envstoch_red, dirpath)Comparing the results of all four simulation runs:
# Get results:
range_df_envstoch_white <- readRange(s_envstoch_white, dirpath)
range_df_envstoch_red <- readRange(s_envstoch_red, dirpath)
par(mfrow=c(2,2),mar=c(3,3,1,1)+.1, tcl=-.1, mgp=c(1.6,.3,0))
# Get common y-range of all 4 simulations
ylim <- c(min(min(range_df$NInds),min(range_df_long$NInds), min(range_df_envstoch_white$NInds), min(range_df_envstoch_red$NInds)), max(max(range_df$NInds),max(range_df_long$NInds), max(range_df_envstoch_white$NInds), max(range_df_envstoch_red$NInds)))
plotAbundance(range_df, ylim = ylim)
plotAbundance(range_df_long, ylim = ylim)
plotAbundance(range_df_envstoch_white, ylim = ylim)
plotAbundance(range_df_envstoch_red, ylim = ylim)
par(mfrow=c(2,2),mar=c(3,3,1,1)+.1, tcl=-.1, mgp=c(1.6,.3,0))
ylim <- c(min(min(range_df$NOccupCells),min(range_df_long$NOccupCells), min(range_df_envstoch_white$NOccupCells), min(range_df_envstoch_red$NOccupCells)), max(max(range_df$NOccupCells),max(range_df_long$NOccupCells), max(range_df_envstoch_white$NOccupCells), max(range_df_envstoch_red$NOccupCells)))
plotOccupancy(range_df, ylim = ylim)
plotOccupancy(range_df_long, ylim = ylim)
plotOccupancy(range_df_envstoch_white, ylim = ylim)
plotOccupancy(range_df_envstoch_red, ylim = ylim)
# Map mean abundance scenario 3
pop_df_envstoch_white <- readPop(s_envstoch_white, dirpath)
ext_envstoch_white <- c(min(pop_df_envstoch_white$x)-500,max(pop_df_envstoch_white$x)+500,min(pop_df_envstoch_white$y)-500,max(pop_df_envstoch_white$y)+500)
r_years_envstoch_white <- stack_pop(pop_df_envstoch_white, ext_envstoch_white, mask=UKmap)
levelplot(r_years_envstoch_white[['mean.NInd.90']], margin=F, scales=list(draw=FALSE), at=c(0,seq(1,max(pop_df_envstoch_white$NInd),length=20)), col.regions=c('grey',rev(inferno(20)))) +
layer(sp.polygons(rasterToPolygons(SpDist, dissolve=F), fill=NA, col='red'))
# Which replicate had highest occupancy?
rep_x = subset(range_df_envstoch_white,Year==95)$Rep[which.max(subset(range_df_envstoch_white,Year==95)$NOccupCells)]
r_years_envstoch_white_1 <- stack_pop(pop_df_envstoch_white, ext_envstoch_white, rep=rep_x, mask=UKmap)
levelplot(r_years_envstoch_white_1[['NInd.90']], margin=F, scales=list(draw=FALSE), at=c(0,seq(1,max(pop_df_envstoch_white$NInd),length=20)), col.regions=c('grey',rev(inferno(20)))) +
layer(sp.polygons(rasterToPolygons(SpDist, dissolve=F), fill=NA, col='red'))
# Map mean abundance scenario 4
pop_df_envstoch_red <- readPop(s_envstoch_red, dirpath)
ext_envstoch_red <- c(min(pop_df_envstoch_red$x)-500,max(pop_df_envstoch_red$x)+500,min(pop_df_envstoch_red$y)-500,max(pop_df_envstoch_red$y)+500)
r_years_envstoch_red <- stack_pop(pop_df_envstoch_red, ext_envstoch_red, mask=UKmap)
levelplot(r_years_envstoch_red[['mean.NInd.90']], margin=F, scales=list(draw=FALSE), at=c(0,seq(1,max(pop_df_envstoch_red$NInd),length=20)), col.regions=c('grey',rev(inferno(20)))) +
layer(sp.polygons(rasterToPolygons(SpDist, dissolve=F), fill=NA, col='red'))
# Which replicate had highest occupancy?
rep_x = subset(range_df_envstoch_red,Year==95)$Rep[which.max(subset(range_df_envstoch_red,Year==95)$NOccupCells)]
r_years_envstoch_red_1 <- stack_pop(pop_df_envstoch_red, ext_envstoch_red, rep=rep_x, mask=UKmap)
levelplot(r_years_envstoch_red_1[['NInd.90']], margin=F, scales=list(draw=FALSE), at=c(0,seq(1,max(pop_df_envstoch_red$NInd),length=20)), col.regions=c('grey',rev(inferno(20)))) +
layer(sp.polygons(rasterToPolygons(SpDist, dissolve=F), fill=NA, col='red'))
2.2 Exercise 2:
Landscape-scale connectivity, matrix permeability and dispersal behaviour
In this second example, RangeShiftR is used at the landscape scale to model functional connectivity of a woodland network for a hypothetical woodland species. The aims are:
- to illustrate how the program can be used to investigate connectivity issues as well as species spatial dynamics at local and landscape scales;
- to show how the program can be run as patch-based;
- to show how additional complexity in the population dynamics and dispersal behaviour can be incorporated;
- and to show how the connectivity analyses can be dependent upon the type of model and on the modelled dispersal behaviour.
We want to reproduce Figure 3 of Bocedi et al. (2014). To this end, we run four different scenarios:
- Explicit sexual model. Constant per-step mortality probability of 0.01. Individuals settle only if at least one individual of the opposite sex is present in the patch (Figure 3b in the paper).
- As in (a), but with different settlement rules. Females settle in suitable patches, while males will settle only if at least one female is present in the patch (Figure 3c in the paper).
- Only-female model. Constant per-step mortality probability of 0.01. Females settle in suitable patches (Figure 3d in the paper).
- As in (a), but with habitat-specific per-step mortality.
As measures for illustrating the connectivity between the initial patch and the rest of the woodland network Fig 3 of Bocedi et al. (2014) uses the ‘final probability of occupancy’ and the ‘mean time to first colonization’. They highlight the dependency of outcomes on the landscape characteristics and movement abilities of the species and, importantly, also on the population dynamics. Notably, these measures both represent multi-generation connectivity.
2.2.1 Create a RS directory
We need to set up the folder structure again with the three sub-folders named ‘Inputs’, ‘Outputs’ and ‘Output_Maps’.
# relative path from working directory:
dirpath = "Tutorials/Exercise2/"
dir.create(paste0(dirpath,"Inputs"), showWarnings = TRUE)
dir.create(paste0(dirpath,"Outputs"), showWarnings = TRUE)
dir.create(paste0(dirpath,"Output_Maps"), showWarnings = TRUE)Copy the map files provided for exercise 2 into the ‘Inputs’ folder.
2.2.2 Landscape parameters
We use a typical British lowland, agricultural landscape having small fragments of woodland, as used by Forest Research, UK, in Watts et al. (2010). The landscape map has an extent of 10km by 6km and a resolution of 10m. Land-covers were aggregated into seven categories (Figure 3a in the paper). Like in the first exercise, the map, landscape_10m_batch.txt, is a raster map with codes for different land-cover types. Land-covers were aggregated into seven categories which are, for simplicity, set as sequential integer numbers:
- 1 = semi-natural broad-leaved woodland
- 2 = planted/felled broad-leaved and mixed woodland, shrubs and bracken
- 3 = heathland, marshy grassland
- 4 = unimproved grassland, mire
- 5 = planted/felled coniferous woodland, semi-improved grassland, swamp
- 6 = improved grasslands, arable, water
- 7 = roads, buildings
landsc <- raster(paste0(dirpath, "Inputs/landscape_10m_batch.txt"))
# Plot land cover map and highlight cells with initial species distribution - option 2 with categorical legend:
landsc.f <- as.factor(landsc)
# add the land cover classes to the raster attribute table
rat <- levels(landsc.f)[[1]]
rat[["landcover"]] <- c("semi-natural broad-leaved woodland", "planted/felled broad-leaved and mixed woodland", "heathland, marshy grassland", "unimproved grassland", "planted/felled coniferous woodland", "improved grasslands, arable, water", "roads, buildings")
levels(landsc.f) <- rat
levelplot(landsc.f, margin=F, scales=list(draw=FALSE), col.regions=brewer.pal(n = 7, name = "Spectral"))
The second text file, woodland_1ha_patchIDs.txt, contains the patch-matrix landscape. It has the same extent and resolution as the land-type map, and each cell contains a unique patch ID that indicates to which patch it belongs. Patch number 0 designates the matrix patch, i.e. all unsuitable habitat.
patch <- raster(paste0(dirpath, "Inputs/woodland_1ha_patchIDs.txt"))
# We can have a glimpse at how many cells the different patches contain:
table(values(patch))##
## 0 1 2 3 4 5 6 7 8 9
## 585734 287 232 243 996 240 238 181 141 990
## 10 11 12 13 14 15 16 17 18 19
## 162 221 311 207 594 694 118 137 172 245
## 20 21 22 23 24 25 26 27 28 29
## 361 423 349 145 1141 138 401 280 336 706
## 30 31 32 33 34 35 36 37 38 39
## 1919 249 154 166 524 215 1277 383 735 113
## 40 41 42 43 44 45 46 47 48 49
## 1008 447 125 100 547 116 225 675 189 110
## 50
## 301# Plot the patches in different colours:
levelplot(patch, margin=F, scales=list(draw=FALSE), at=0:50, colorkey=F,
col.regions = c('black',rep(brewer.pal(n = 12, name = "Paired"),5))
) 
The last text file, patch30.txt, is a map that specifies the patch(es) that contain(s) the initial distribution of the species. In our case, this is only the patch with ID 30.
patch30 <- raster(paste0(dirpath, "Inputs/patch30.txt"))
# Look at initial patch:
plot(patch30)
We are ready to set up the landscape parameter object with these maps, their respective resolutions, and the carriying capacities for all land types. For the latter we choose to define only ‘semi-natural broad-leaved woodland’ (code 1) as suitable for our species.
land <- ImportedLandscape(LandscapeFile = "landscape_10m_batch.txt",
PatchFile = "woodland_1ha_patchIDs.txt",
Resolution = 10,
Nhabitats = 7,
K = c(10, rep(0,6)),
SpDistFile = "patch30.txt",
SpDistResolution = 10)2.2.3 Scenario a) Demographic and dispersal parameters
Sexual model with both sexes settling only in the presence of a mate
We will simulate a sexual species with simple, two-staged stage-structured population dynamics. The parameters are chosen to be representative of species having moderately high fecundity, high juvenile mortality and low adult mortality. This is encoded in the following transition matrix
## [,1] [,2] [,3]
## [1,] 0 0.0 5.0
## [2,] 1 0.1 0.0
## [3,] 0 0.4 0.8The first row and column describe the juvenile stage, the others the two adult stages. Juveniles will develop to the first adult stage at the end of their first year with a probability of 1.0, which allows for juvenile dispersal before any mortality happens.
In order to add a stage-structure to our population dynamics, we use the StageStructure function within the demography module.
The reproduction type 1 denotes a simple sexual model.
trans_mat <- matrix(c(0, 1, 0, 0, 0.1, 0.4, 5, 0, 0.8), nrow = 3, byrow = F)
demo <- Demography(ReproductionType = 1, # simple sexual model
StageStruct = StageStructure(Stages=3, # 1 juvenile + 2 adult stages
TransMatrix=trans_mat,
MaxAge=1000,
SurvSched=2,
FecDensDep=T))After reproduction, only juveniles shall be able to disperse, and this according to a density-dependent emigration probability. Therefore we enable the options DensDep and StageDep and in the matrix EmigProb we set the parameters D0 = 0.5, α = 10.0 and β = 1.0 for juveniles and all to zero for the other (i.e. all adult) stages.
To account for functional connectivity, we use a mechanistic movement model which enables individuals to interact with the landscape and determine their path according to what they can perceive in the landscape. Therefore we will simulate movements with a stochastic movement simulator (SMS) where individuals move stepwise (each step being one cell) and the direction chosen at each step is determined by the land cover costs (specified for each land type), the species’ perceptual range (PR) and directional persistence (DP). We set these parameters so that individuals have a perceptual range of 50m, use the arithmetic mean method (the default) for calculating effective cost (which tends to emphasize the avoidance of high-cost landscape features) and tend to follow highly correlated paths within the landscape. We also set a constant per-step mortality probability (StepMort).
Once arrived in a new patch, an individual can decide to settle or not based on certain settlement rules. Finding suitable habitat is in any case a necessary condition. Additionally, we set the mating requirement, i.e. there has to be at least one individual of the opposite sex present for the patch to be considered suitable for settlement.
disp <- Dispersal(Emigration = Emigration(DensDep=T, StageDep=T,
EmigProb = cbind(0:2,c(0.5,0,0),c(10.0,0,0),c(1.0,0,0)) ),
Transfer = SMS(PR=5, DP=10, Costs = c(1,1,3,5,10,20,50), StepMort = 0.01),
Settlement = Settlement(FindMate = T) )2.2.4 Scenario a) Initialisation & simulation
We choose to initialise our simulation in all the patches that are specified by the loaded species distribution map (in our case only patch number 30) at a density of 10 individuals per hectare with an equal number of individuals in stages 1 and 2 and with their respective minimum age.
# Population is initialised in Patch 30:
init <- Initialise(InitType = 1, # from loaded species distribution map
SpType = 0, # all suitable cells
InitDens = 2, # user-specified density
IndsHaCell = 10,
PropStages = c(0,0.5,0.5),
InitAge = 0)We set the simulation time to 100 years and 20 replicates, and set the output types to write the files for population, occupancy and range data every year.
sim <- Simulation(Simulation = 0,
Replicates = 20,
Years = 100,
OutIntPop = 1,
OutIntOcc = 1,
OutIntRange = 1)As before, we need to put all modules together to a parameter master.
s <- RSsim(batchnum = 3, land = land, demog = demo, dispersal = disp, simul = sim, init = init)And run the simulation:
RunRS(s, dirpath)2.2.5 Scenario a) Analyse output
To analyse the simulation output, we first plot the meta-population results. Note here that - in contrast to the cell-based model from exercise 1 - the plotted occupancy refers to occupied patches rather than cells.
par(mfrow=c(1,2))
plotAbundance(s, dirpath)
plotOccupancy(s, dirpath)
In order to create occupancy maps, we first plot the landscape with the suitable patches in green and the initial patch in red. This color scheme was also used in Fig. 3a of Bocedi et al. (2014).
# We have initiated the population in the patch with ID=30. We highlight this in the map.
values(patch30)[values(patch30)<1] <- NA
values(patch)[values(patch)<1] <- NA
levelplot(landsc, margin=F, scales=list(draw=FALSE),at=seq(.5,7.5,by=1), colorkey=F,
col.regions = rev(brewer.pal(n = 7, name = "Greys") )) +
levelplot(patch, margin=F, scales=list(draw=FALSE), col.regions="green4") +
layer(sp.polygons(rasterToPolygons(patch30, dissolve=T), fill=NA, col='red',lwd=2))
# Store underlying landscape map display for later:
bg <- function(main=NULL){
levelplot(landsc, margin=F, scales=list(draw=FALSE),at=seq(.5,7.5,by=1), colorkey=F,
col.regions = rev(brewer.pal(n = 7, name = "Greys") ), main=main)
}To reproduce Fig. 3b of Bocedi et al. (2014), we map the mean occupancy probability for each patch after 100 years (left panel in the paper) as well as the mean time to colonisation (right panel), both calculated over the 20 replicates. We need to do a little data processing for this.
# read population results
pop_df <- readPop(s, dirpath)
pop_df <- pop_df[pop_df$PatchID!=0,]
# Occupancy probability after 100 years over all replicates
occ_prob <- sapply(sort(unique(pop_df$PatchID)), FUN=function(x,n=length(unique(pop_df$Rep))){ind <- sum(subset(pop_df,PatchID==x & Year==100)$NInd>0)/n; names(ind)=x;ind})
# Time to colonisation:
col_time <- data.frame(sapply(sort(unique(pop_df$PatchID)),FUN=function(p){sapply(sort(unique(pop_df$Rep)),FUN=function(r){ ifelse(nrow(subset(pop_df,PatchID==p & Rep==r & NInd>0)), min(subset(pop_df,PatchID==p & Rep==r & NInd>0)$Year), NA) })}))
names(col_time) <- sort(unique(pop_df$PatchID))
col_time_mean <- colMeans(col_time)
# Update patch information
patch_occ_prob <- patch_col_time <- patch
values(patch_occ_prob)[values(patch)>0] <- 0
values(patch_col_time)[values(patch)>0] <- -9
for (i in as.numeric(names(occ_prob))){
values(patch_occ_prob)[values(patch)==i] <- occ_prob[paste(i)]
}
for (i in as.numeric(names(col_time_mean))){
values(patch_col_time)[values(patch)==i] <- ifelse(is.na(col_time_mean[paste(i)]),-9,col_time_mean[paste(i)])
}
#values(patch_occ_prob)[values(patch)<1] <- values(patch_col_time)[values(patch)<1] <- NA
# map occupancy probability
mycol_occprob <- colorRampPalette(c('blue','orangered','gold'))
levelplot(patch_occ_prob, margin=F, scales=list(draw=FALSE), at=seq(0,1,length=11), col.regions=mycol_occprob(11))
# map occupancy probability on landscape background. For this, we first define a colorkey function
col.key <- function(mycol, at, space='bottom',pos=0.05, height=0.6, width=1) {
key <- draw.colorkey(
list(space=space, at=at, height=height, width=width,
col=mycol)
)
key$framevp$y <- unit(pos, "npc")
return(key)
}
bg() + levelplot(patch_occ_prob, margin=F, scales=list(draw=FALSE), at=seq(0,1,length=11), col.regions=mycol_occprob(11))
grid.draw(col.key(mycol_occprob(11),at=seq(0,1,length=11)))
# map colonisation time
mycol_coltime <- colorRampPalette(c('orangered','gold','yellow','PowderBlue','LightSeaGreen'))
levelplot(patch_col_time, margin=F, scales=list(draw=FALSE), at=c(-9,seq(-.001,max(pop_df$Year),length=11)), col.regions=c('blue',mycol_coltime(11)))
# map colonisation time on landscape background
bg() + levelplot(patch_col_time, margin=F, scales=list(draw=FALSE), at=c(-9,seq(-.001,max(pop_df$Year),length=11)), col.regions=c('blue',mycol_coltime(11)))
grid.draw(col.key(c('blue',mycol_coltime(11)), c(-9,seq(-.001,max(pop_df$Year),length=11))))
For convenience, let’s put all the population data processing into a function for later usage.
# Put all of this into a function
get_patch_results <- function(pop_df,patch_r){
# This function takes the patch-based population file and the patch raster map and produces summary stats and maps on occupancy probability and mean time to colonisation
# Occupancy probability after 100 years over all replicates
occ_prob <- sapply(sort(unique(pop_df$PatchID)), FUN=function(x,n=length(unique(pop_df$Rep))){ind <- sum(subset(pop_df,PatchID==x & Year==100)$NInd>0)/n; names(ind)=x;ind})
# Time to colonisation:
col_time <- data.frame(sapply(sort(unique(pop_df$PatchID)),FUN=function(p){sapply(sort(unique(pop_df$Rep)),FUN=function(r){ ifelse(nrow(subset(pop_df,PatchID==p & Rep==r & NInd>0)), min(subset(pop_df,PatchID==p & Rep==r & NInd>0)$Year), NA) })}))
names(col_time) <- sort(unique(pop_df$PatchID))
col_time_mean <- colMeans(col_time)
patch_occ_prob <- patch_col_time <- patch
values(patch_occ_prob)[values(patch)>0] <- 0
values(patch_col_time)[values(patch)>0] <- -9
for (i in as.numeric(names(occ_prob))){
values(patch_occ_prob)[values(patch)==i] <- occ_prob[paste(i)]
}
for (i in as.numeric(names(col_time_mean))){
values(patch_col_time)[values(patch)==i] <- ifelse(is.na(col_time_mean[paste(i)]),-9,col_time_mean[paste(i)])
}
values(patch_occ_prob)[values(patch)<1] <- values(patch_col_time)[values(patch)<1] <- NA
return(list(occ_prob=occ_prob, col_time=col_time, col_time_mean=col_time_mean, patch_occ_prob=patch_occ_prob, patch_col_time=patch_col_time))
}
summary_pop_a <- get_patch_results(pop_df,patch)2.2.6 Scenario b)
Females settle independent of males
This experiment was designed to provide an example of how the dispersal behaviour of the species and the specification of settlement rules can change the estimated connectivity of a habitat network. We will relax the mating requirement a little by making it sex-dependent and only activating it for males. This means that female dispersers will settle in suitable patches regardless of males, while males settle only when finding a mate.
# Change Settlement rules
disp_b <- Dispersal(Emigration = Emigration(DensDep=T, StageDep=T,
EmigProb = cbind(0:2,c(0.5,0,0),c(10.0,0,0),c(1.0,0,0)) ),
Transfer = SMS(PR=5, DP=10, Costs = c(1,1,3,5,10,20,50), StepMort = 0.01),
Settlement = Settlement(FindMate = c(F,T), SexDep=T, Settle=cbind(c(0,1)) ) )
# Update simulation
sim_b <- Simulation(Simulation = 1,
Replicates = 20,
Years = 100,
OutIntPop = 1,
OutIntOcc = 1,
OutIntRange = 1)
# Update parameter master
s_b <- s + disp_b + sim_b# run simulation
RunRS(s_b, dirpath)Now, let’s post-process the simulation results and plot the maps.
# Read in population results
pop_df_b <- readPop(s_b, dirpath)
# Process the results to get occupancy probabilities and colonisation times
summary_pop_b <- get_patch_results(pop_df_b,patch)
# Map occupancy probabilities:
bg() + levelplot(summary_pop_b$patch_occ_prob, margin=F, scales=list(draw=FALSE), at=seq(0,1,length=11), col.regions=mycol_occprob(11))
grid.draw(col.key(mycol_occprob(11),at=seq(0,1,length=11)))
# map colonisation time + background
bg() + levelplot(summary_pop_b$patch_col_time, margin=F, scales=list(draw=FALSE), at=c(-9,seq(-.001,max(pop_df$Year),length=11)), col.regions=c('blue',mycol_coltime(11)))
grid.draw(col.key(c('blue',mycol_coltime(11)), c(-9,seq(-.001,max(pop_df$Year),length=11))))
From both the visualisation and the results, you will be able to see how this change in the settlement rules substantially increases the number of patches which in at least some simulations become occupied, their probability of occupancy and the mean time to colonization, changing the overall modelled functional connectivity of the woodland network over 100 years.
2.2.7 Scenario c)
Asexual / Female-only model
Here we change the way we model the population dynamics such that we use an only-female model. This change also has important consequences for the dispersal process and potential implications for patterns of colonization across a landscape. The stage-structured model remains the same apart from being an only-female model. This assumes that males are not limiting, and that the population dynamics are driven only by females. It also means that sexes are not modelled explicitly and it is not possible to account for behaviours like mate-finding in the settlement decisions; females will settle in suitable habitat patches and then will automatically be able to attempt reproduction.
Because now only females are considered, some parameters need to be changed. We set the fecundity of stage 3 to 2.5 instead of 5.0 and 1/b to 5 instead of 10. Sex-dependent settlement options are no longer available.
# Change carrying capacitiy to half its value of the sexual model
land_c <- ImportedLandscape(LandscapeFile = "landscape_10m_batch.txt",
PatchFile = "woodland_1ha_patchIDs.txt",
Resolution = 10,
Nhabitats = 7,
K = c(5, rep(0,6)),
SpDistFile = "patch30.txt",
SpDistResolution = 10)
# Change demography settings
demo_c <- Demography(ReproductionType = 0, # female-only model
StageStruct = StageStructure(Stages=3,
TransMatrix=matrix(c(0, 1, 0, 0, 0.1, 0.4, 2.5, 0, 0.8), nrow = 3, byrow = F),
MaxAge=1000,
SurvSched=2,
FecDensDep=T))
# Remove settlement rules
disp_c <- Dispersal(Emigration = Emigration(DensDep=T, StageDep=T,
EmigProb = cbind(0:2,c(0.5,0,0),c(10.0,0,0),c(1.0,0,0)) ),
Transfer = SMS(PR=5, DP=10, Costs = c(1,1,3,5,10,20,50), StepMort = 0.01),
Settlement = Settlement()
)
# Update simulation
sim_c <- Simulation(Simulation = 2,
Replicates = 20,
Years = 100,
OutIntPop = 1,
OutIntOcc = 1,
OutIntRange = 1)
# parameter master
s_c <- RSsim(batchnum = 3, land = land_c, demog = demo_c, dispersal = disp_c, simul = sim_c, init = init)RunRS(s_c, dirpath)We do the same output processing as before and plot the occupancy maps.
# Read population results
pop_df_c <- readPop(s_c, dirpath)
summary_pop_c <- get_patch_results(pop_df_c,patch)
# Map occupancy probabilities:
bg() + levelplot(summary_pop_c$patch_occ_prob, margin=F, scales=list(draw=FALSE), at=seq(0,1,length=11), col.regions=mycol_occprob(11))
grid.draw(col.key(mycol_occprob(11),at=seq(0,1,length=11)))
# Map colonisation time + background
bg() + levelplot(summary_pop_c$patch_col_time, margin=F, scales=list(draw=FALSE), at=c(-9,seq(-.001,max(pop_df$Year),length=11)), col.regions=c('blue',mycol_coltime(11)))
grid.draw(col.key(c('blue',mycol_coltime(11)), c(-9,seq(-.001,max(pop_df$Year),length=11))))
As we see from the results, not accounting explicitly for sexes and settlement behaviours leads to a drastic increase in the overall occupancy of the habitat network after 100 years.
2.2.8 Scenario d)
Habitat-specific per-step mortality
In this last experiment, we will demonstrate how RangeShiftR can incorporate more complexity in the way that movement is modelled. We relax the unrealistic assumption that the per-step mortality is constant across all the land-cover types, and assign different mortality values to each habitat. To set up this simulation, we use the parameters from scenario a) and only add a modified transfer module. Here, set StepMort habitat-dependent by giving it a vector that specifies the mortality probabilities for each land-type.
# Update Transfer sub-module within the dispersal module
disp_d <- disp + SMS(PR=5, DP=10, Costs = c(1,1,3,5,10,20,50),
StepMort = c(0,0,0,0.01,0.01,0.02,0.05)
)
# Update simulation
sim_d <- Simulation(Simulation = 3,
Replicates = 20,
Years = 100,
OutIntPop = 1,
OutIntOcc = 1,
OutIntRange = 1)
# Use parameter master from a) and add new transfer module
s_d <- s + disp_d + sim_dRunRS(s_d, dirpath)As before, we proceed with the known output processing.
# Read population results
pop_df_d <- readPop(s_d, dirpath)
summary_pop_d <- get_patch_results(pop_df_d,patch)
# Map occupancy probabilities:
bg() + levelplot(summary_pop_d$patch_occ_prob, margin=F, scales=list(draw=FALSE), at=seq(0,1,length=11), col.regions=mycol_occprob(11))
grid.draw(col.key(mycol_occprob(11),at=seq(0,1,length=11)))
# map colonisation time + background
bg() + levelplot(summary_pop_d$patch_col_time, margin=F, scales=list(draw=FALSE), at=c(-9,seq(-.001,max(pop_df$Year),length=11)), col.regions=c('blue',mycol_coltime(11)))
grid.draw(col.key(c('blue',mycol_coltime(11)), c(-9,seq(-.001,max(pop_df$Year),length=11))))
We see that such small changes in the per-step mortality, in interaction with the landscape structure, make a big difference in the results, in this case decreasing the functional connectivity of the network.
2.2.9 Summary
Let’s plot all maps next to each other.
# Plot occupancy probabilities for all scenarios
occ_a <- bg(main="Scenario A") + levelplot(summary_pop_a$patch_occ_prob, margin=F, scales=list(draw=FALSE), at=seq(0,1,length=11), col.regions=mycol_occprob(11))
occ_b <- bg(main="Scenario B") + levelplot(summary_pop_b$patch_occ_prob, margin=F, scales=list(draw=FALSE), at=seq(0,1,length=11), col.regions=mycol_occprob(11))
occ_c <- bg(main="Scenario C") + levelplot(summary_pop_c$patch_occ_prob, margin=F, scales=list(draw=FALSE), at=seq(0,1,length=11), col.regions=mycol_occprob(11))
occ_d <- bg(main="Scenario D") + levelplot(summary_pop_d$patch_occ_prob, margin=F, scales=list(draw=FALSE), at=seq(0,1,length=11), col.regions=mycol_occprob(11))
grid.arrange(occ_a,occ_b,occ_c,occ_d, ncol=2)
grid.draw(col.key(mycol_occprob(11),at=seq(0,1,length=11), space='right',pos=0.5))
# Plot colonisation times for all scenarios
# colonisation times
col_a <- bg(main="Scenario A") + levelplot(summary_pop_a$patch_col_time, margin=F, scales=list(draw=FALSE), at=c(-9,seq(-.001,max(pop_df$Year),length=11)), col.regions=c('blue',mycol_coltime(11)))
col_b <- bg(main="Scenario B") + levelplot(summary_pop_b$patch_col_time, margin=F, scales=list(draw=FALSE), at=c(-9,seq(-.001,max(pop_df$Year),length=11)), col.regions=c('blue',mycol_coltime(11)))
col_c <- bg(main="Scenario C") + levelplot(summary_pop_c$patch_col_time, margin=F, scales=list(draw=FALSE), at=c(-9,seq(-.001,max(pop_df$Year),length=11)), col.regions=c('blue',mycol_coltime(11)))
col_d <- bg(main="Scenario D") + levelplot(summary_pop_d$patch_col_time, margin=F, scales=list(draw=FALSE), at=c(-9,seq(-.001,max(pop_df$Year),length=11)), col.regions=c('blue',mycol_coltime(11)))
grid.arrange(col_a,col_b,col_c,col_d, ncol=2)
grid.draw(col.key(c('blue',mycol_coltime(11)), c(-9,seq(-.001,max(pop_df$Year),length=11)), space='right',pos=0.5))
References
Bocedi, G., S.C.F. Palmer, G. Pe’er, R.K. Heikkinen, Y.G. Matsinos, K. Watts, and J.M.J. Travis. 2014. “RangeShifter: A Platform for Modelling Spatial Eco-Evolutionary Dynamics and Species’ Responses to Environmental Changes.” Methods in Ecology and Evolution 5 (4): 388–96. https://doi.org/10.1111/2041-210X.12162.
Morton, D., C. Rowland, C. Wood, L. Meek, C. Marston, G. Smith, R. Wadsworth, and I.. Simpson. 2011. “Final Report for Lcm2007 - the New Uk Land Cover Map.” Countryside Survey Technical Report No 11/07. NERC/Centre for Ecology & Hydrology.